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PHASE  CHANGE  AROUND  A  CIRCULAR  PIPE 

V.J.  Lunardini 


INTRODUCTION 

The  question  of  freezing  and  thawing  of  soils  is  im¬ 
portant  for  engineering  design  in  permafrost  regions. 

For  example,  the  support  of  buried  pipes,  utilidors, 
and  foundations  often  relies  on  the  mechanical 
strength  of  the  permafrost  and  it  is  necessary  to  know 
the  extent  of  thaw  endured  in  the  permafrost  by  heated 
structures.  Thus  solutions  are  sought  to  conduction 
problems  with  moving  boundaries  for  analysis  of  these 
questions  and  synthesis  of  solutions.  Interest  in  these 
types  of  problems  has  also  been  stimulated  by  the  re¬ 
quirements  of  latent  heat  storage  in  solar  energy 
systems. 

For  a  homogeneous,  semi-infinite  slab,  initially  at  a 
temperature  different  than  the  fusion  value,  there 
exists  a  general  solution  for  a  step  change  of  surface 
temperature  (Neumann  1860).  Recently,  Lunardini 
and  Varotta  (1980)  have  given  an  accurate,  but  approx¬ 
imate,  solution  to  this  problem  in  convenient  form  by 
using  the  heat  balance  integral  technique.  This  solu¬ 
tion  is  reproduced  here  since  it  will  be  used  later  for 
the  cylindrical  problem.  The  phase  change  depth,  for 
the  semi-infinite  medium,  is  given  by 

X  =  2  yy/^t  (1) 

72  _  ^ 

2  a 

where 

a=  (sT+2+^£t)(ST+2) 

V  a21  / 

6,  =  -25-,  (Sj ♦ 2 ♦  k*L  0il )  -1 1*11? 5tI2. 

\  «2,  /  3  a21 

ST  =c,(rj-Tf)/8 


<t>=(T,-T0)l(Ts-Tf). 

The  superheat  parameter  <t>  is  a  measure  of  the  amount 
the  initial  temperature  differs  from  the  fusion  temper¬ 
ature.  ST  and  <t>  are  for  the  case  of  thawing.  Freezing 
relations  are  givgn  in  the  nomenclature.  The  Stefan 
number  Sy  is  the  ratio  of  the  sensible  to  latent  heat, 
while  <t>  is  the  superheat  or  subcooling  parameter. 

In  contrast  to  the  Neumann  problem,  there  are  no 
general,  exact  solutions  for  phase  change  in  cylindrical 
coordinates.  Since  this  is  an  important  geometry  for 
practical  applications,  a  significant  effort  has  been  ex¬ 
pended  upon  analytical  solutions  limited  to  certain  do¬ 
mains  or  approximate  solutions.  The  geometry  of  the 
system  (Fig.  1)  is  that  of  a  cylinder  surrounded  by  an 
infinite,  homogeneous  medium.  A  review  will  be  made 
of  several  simple  cases  followed  by  new  solutions  of 
the  general  problem. 


ZERO  SUPERHEAT,  0  =  0 


A  number  of  solutions  are  available  for  the  case 
when  the  temperature  of  the  medium  is  initially  uni¬ 
form  at  the  fusion  temperatures.  This  condition  greatly 
simplifies  the  problem  since  only  one  phase  need  be 
considered;  that  which  is  changing  phase. 

The  energy  condition  at  the  phase  change  interface 
is  given  by  the  same  relation  as  in  the  cartesian  case 


-k^h +k2dh  =  ±P^A 

'dr  2  dr  dt 


(3) 


where  the  upper  sign  is  for  melting  and  the  lower  is  for 
freezing.  This  assumes  that  the  interface  motion  is  in 
the  positive  direction  of  the  coordinate  r,  and  that  the 
density  is  constant. 


Figure  1.  Geometry  of  cylindrical 
system. 


Constant  phase  change  rate 

A  simple  system  which  has  an  exact  solution  re¬ 
quires  that  the  temperature  gradient  at  the  phase  change 
interface  be  a  constant.  This  case  has  been  discussed  by 
Kreith  and  Romie  (1955)  for  inward  solidification.  The 
problem  for  outward  solidification  can  be  written  as 
follows: 


LA/z-iTUL^I 

(4) 

r  3r\  dr)  a  3 1 

T(R,  t)  =  rf 

(4a) 

(3=  U 


G'o 

(Tr  rs) 


r 


(6) 


where 


T  =  al  ! 


Kreith  and  Romie  (1955)  used  an  iterative-type  series 
solution  for  the  temperature  of  the  solid.  This  will  not 
be  given  here  since  the  problem  is  not  too  practical  due 
to  the  necessity  of  imposing  a  complicated  transient 
temperature  at  the  cylinder  surface. 

Zero  sensible  heat,  5T  =  0 

If  the  sensible  heat  of  the  material  is  also  neglected, 
then  particularly  simple  solutions  are  available.  Since 
the  Stefan  number  is  a  measure  of  the  sensible  heat, 
this  situation  is  equivalent  to  assuming  that  the  Stefan 
number  is  zero.  Carslawand  Jaeger  (1959)  presented 
a  thaw  solution  for  the  case  using  the  quasi-steady  ap¬ 
proximation  when  the  surrounding  medium  is  at  the 
phase  change  temperature  Tf.  The  quasi-steady  method 
is  accurate  if  5T  «  1 .0.  The  problem  is  as  follows: 


kmA-P^R 

dr  dt 

(4b) 

o 

II 

Q-. _ 

(7) 

T(ro,0)  =  Tf. 

(4c) 

T=  Tf;  r  =  R 

(7a) 

The  temperature  gradient  at  the  phase  change  inter- 

r=r5;  r  =  rQ. 

(7b) 

face  is  specified  as 


Equation  3,  the  interface  balance,  is  now 


It  is  not  necessary  to  solve  eq  4  in  order  to  find  the 

location  of  the  freezing  front.  From  eq  4b  and  4d  The  solution  for  the  phase  change  location  is  straight¬ 

forward  and  is 


dR  _  k  r 

Vt'Tt 


(5) 


ip2  lnp-@—  +1=  r. 
2  4  4 


(8) 


Thus  the  rate  of  movement  of  the  phase  change  inter¬ 
face  is  a  constant.  The  location  is 


R  = 


*Jf+r 

Pt 


O 


Equation  8  will  overestimate  the  thaw  beneath  the  pipe 
since  all  of  the  energy  transfer  from  the  pipe  will  go 
into  thawing.  A  correction  can  be  made  using  an  effec¬ 
tive  latent  heat  in  place  of  8.  One  example  of  an  effec¬ 
tive  latent  heat  is 


8e  =  «(l  +  C2,0ST  +  ^). 


2 


The  surface  heat  transfer  is 


The  heat  transfer  at  the  cylinder  surface  can  then  be 
written 

211/8,6  ln/3 

Finite  sensible  heat 

The  heat  balance  integral  can  also  be  conveniently 
used  for  the  no  superheat  problem  and  may  also  in¬ 
clude  large  values  of  the  Stefan  number.  The  equa- 


tions  are 

rot*0 

ll 

(9) 

T(R,  t)  =  Tf 

(9a) 

T{rv  t)  =  rs 

(9b) 

H 

©■ 

(9c) 

-M*T&‘)=ptdR  , 
dr  dt 

(9d) 

Integration  of  eq  9  once,  over  the  space  coordinate, 
yields  the  heat  balance  integral  equation 


[  dr  0  dr 


=  M_RTdR 

dt  dt 


(10) 


where  the  integrated  temperature  is 


e  = 


rTdr. 


(ID 


The  solution  method  consists  of  guessing  an  approxi¬ 
mation  for  the  temperature  which  will  satisfy  all  of 
the  conditions  in  terms  of  an  unknown  parameter 
such  as  R  (t).  Lardner  and  Pohle  (T961 )  noted  that  a 
logarithmic  temperature  approximation  is  more  appro¬ 
priate  than  a  polynomial  in  r  since  the  area  is  varying 
with  r.  They  suggest  that  T  =  f[r)  Inr  be  used  as  an 
approximation.  The  temperature  is  thus  assumed  to 
be  a  logarithmic  relation,  satisfying  eq  9a  and  9b. 


T  = 


V 


(Ts-Tf) 


Hr/r0) 
In  (R/r0)' 


(12) 


Equations  10-12  then  yield  a  differential  equation  for 
the  phase  change  location 

°3) 

Equation  1 3  can  be  integrated  to  give 

M^-l  -2(ln0)-(ln0)2-...  + 

4  l  nn!  J 

1  ln0-£+l  =r  (14) 

2  4', 

This  equation  reduces  to  the  zero  sensible  heat  solution 
(eq  8)  when  the  Stefan  number  is  zero.  Equation  14  is 
comparable  in  accuracy  to  numerical  solutions,  valid 
for  0  =  0,  with  any  value  of  ST. 

FINITE  SUPERHEAT 

The  more  practical  problems,  in  which  the  initial 
temperature  is  not  at  the  fusion  value,  present  signifi¬ 
cantly  more  difficult  analyses.  Several  methods  have 
been  utilized.  All  of  the  problems  assume  that  the  pipe 
is  buried  at  an  infinite  depth  because  a  finite  burial 
depth  will  effect  very  severe  restraints  on  the  solution 
methods. 

Quasi-steady  solution 

A  simple  solution  can  again  be  obtained  with  the 
quasi-steady  state  approximation.  Khakimov  (1957) 
investigated  this  problem  and  introduced  the  concept 
of  a  thermal  layer  of  influence  around  the  buried  pipe. 
Consider  the  case  of  thawing  of  a  medium  initially  fro¬ 
zen  atT=T0.  A  hot  cylinder  with  a  surface  tempera¬ 
ture  rs  is  inserted  in  the  medium  at  time  zero.  The 
thaw  effect  (temperature  change)  is  assumed  to  extend, 
at  any  time  t,  to  a  finite  distance  6.  That  is,  the  tem¬ 
perature  of  the  medium  will  be  T0  at  some  location  suf¬ 
ficiently  far  from  the  hot  cylinder  (see  Fig.  1).  This 
concept  of  a  temperature  penetration  is  also  used  in 
the  heat  balance  integral  method  and  in  boundary 
layer  theory. 

From  experimental  evidence,  Khakimov  (1957)  as¬ 
sumed  that  the  ratio  b  =  ( 5/R )  was  a  constant  equal  to 
4.5.  The  assumption  that  b  is  a  constant,  for  a  given 
5t  and  0,  is  correct  for  the  Neumann  problem  (lunar  - 
dini  and  Varotta  1980),  but  is  invalid  for  a  cylindrical 
system.  Nevertheless  it  does  simplify  the  equations 
and  yields  reasonable  approximations. 

The  temperature  will  be  the  solution  of 


3 


±[rdJL\- 0 
dr  \  dr) 


=  dJl 

dt 


(20) 


05.). 


q  =  -k 


o 


for  each  region  with  the  boundary  conditions 


Combining  eq  19  and  20  yields 


«,  =  w2  =  u  f ;  r=R 

(15a) 

-P-  [dt~  /  mx  In*  +x  -  ^  x2-1  dx 
r*  C]  JQ  2  4xlnz 

u\  =us>  r  =  ro 

(15b) 

u2  =  0;  r  =  6 

(15c) 

where 

where 

L-p  +  u£y 
m  - - 

u  =  T-  Ta. 

c,(Ts-rf) 

The  solutions  for  the  temperatures  are 

The  solution  to  this  equation  is 

(7>-rs) 

=  M«W  ***'"“> 

(16) 

4t  =  2mST 02  ln/3  +ST (1  -m)(/12-1)-ST 

“2  = 


uf 

In  (R/8) 


I  n  (/■/  5 )  . 


(17) 


The  total  energy  added  to  the  thawing  medium  will 
be  the  latent  heat  needed  to  thaw  the  layer  between 
r0  and  R  and  the  sensible  heat  used  to  increase  the 
temperature  of  all  the  layers  up  to  r  =  6.  Thus, 


where 


x  In* 


Q  =  t,(R2-^)L+2tiC] 


ru ,  dr  + 


2nC2 


(18) 


Carrying  out  the  integration  yields 


where  L  =  volumetric  latent  heat 
C,,  C2  :  volumetric  specific  heat  of  thawed  and 


frozen  layers 


During  a  small  time  increment  the  change  in  the  energy 
absorbed  by  the  system  must  equal  the  energy  added 
to  the  system  at  the  cylinder  surface. 


The  universal  function  0  has  been  numerically  evalu¬ 
ated  and  is  plotted  as  Figure  2,  allowing  the  thaw  for 
any  soil  to  be  approximated.  The  method  is  simple 
but  limited  in  application  to  first  estimates  of  thaw 
depth.  When  0  is  small  the  method  underpredicts  the 
phase  change  depth  and  the  accuracy  decreases  with 
time  since  b  is  not  constant. 

Heat  balance  integral  solution 

The  heat  balance  integral  method  may  be  applied 
to  the  problem  of  finite  superheat  but  the  labor  is  in¬ 
creased  considerably.  The  following  analysis  will  allow 
simple  numerical  techniques  to  be  applied  to  systems 
with  typical  soil  parameters.  The  problem  may  be 
stated  for  a  melting  system  as  follows: 


d_(r  £Il\  =  r  3ri 
dr\  9r  /  a,  d t 

(23) 

T\  (R,  t)  =  rf 

(23a) 

r 

M 

(23b) 

T,(r,0)  =  T0 

(23c) 

4 


0 

Figure  2.  Quasi-steady  phase  change  parameter. 
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The  heat  balance  integral  equations  are 


Equations  25,  26,  29,  and  30  applied  to  eq  27  and  28 
now  lead  to  the  following  coupled  system  of  equations. 


a125T—  =  — L  ♦  — L 
1  dr  12-1  Infl 

=J_. 

I  U  In*  40(ln/3)2/  J  *  '"0 


R  *h(R,t) 

32"]  (r0,  t 
~r0  ■  -  . 

**<(- 

dr 

dr 

dt  '  dt 
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Table  1 .  Frozen  and  thawed  properties  ratios  used 
with  Figures  6  through  22. 


Thaw  Freeze 


XS 

*12  =  V*f 

ai2  =  <*,/af 

*12  =  *f/*t 

<*12  ~  <*f /« 

0 

1 

1 

1 

1 

0.1 

0.873 

0.792 

1.1455 

1.2626 

0.2 

0.7621 

0.6326 

1.3122 

1.5808 

0.3 

0  6653 

0.5091 

1.5031 

1.9642 

0.4 

0.5808 

0.4121 

1.7218 

2.4266 

0.5 

0.5070 

0.3354 

1.9724 

2.9815 

0.6 

0.4426 

0.2742 

2.2593 

3.6470 

0.7 

0.3864 

0.2252 

2.5880 

4.4405 

0.8 

0.3373 

0.1855 

2.9647 

5.3908 

0.9 

0.2945 

0.1533 

3.3956 

6.5231 

1.0 

0.2571 

0.1271 

3.8895 

7.8685 

where  Cl  =  8/R.  Equations  31  -33  were  solved  numeri¬ 
cally  with  a  fourth  order,  Runge-Kutta,  predictor- 
corrector  technique.  Since  the  problem,  as  specified, 
is  initially  singular  at  the  origin,  the  Neumann  solution 
was  used  to  start  the  calculation.  Sparrow  et  al.  (1978) 
solved  this  problem  numerically,  with  a12  =  /?t2  =  1 


for  a  range  of  ST  and  <t>,  also  using  the  Neumann  solu¬ 
tion  as  a  start.  The  results  of  the  much  simpler  method 
presented  here  are  within  5%  of  the  values  reported  by 
Sparrow  et  al.  The  comparisons  are  shown  in  Figures 
3-5.  The  calculations  have  been  generalized  for  a  range 
of  a, 2,  k  )2  whi  ;h  are  pertinent  for  soil  systems,  and 
are  presented  in  Figures  6-22.  The  property  values 
associated  with  each  value  ofx6  are  given  in  Table  1. 
These  property  ratios  are  obtained  using  the  method  of 
Lunardini  and  Varotta  (1980).  While  the  curves  have 
been  developed  particularly  for  soils,  they  are  valid  for 
any  medium  with  the  same  properties. 

Tien  and  Churchill  (1965)  also  numerically  evalu¬ 
ated  freezing  outside  a  cylinder.  Their  calculations  are 
more  extensive  than  Sparrow  et  al.  (1978)  and  the 
numerical  technique  was  entirely  different  but  the  re¬ 
sults  are  essentially  the  same. 

Approximate  methods 

The  numerical  solutions,  while  very  good,  are  often 
not  as  convenient  as  analytical  solutions.  Thus  further 
work  will  be  carried  out  to  obtain  approximate  solutions 
which  yield  results  with  acceptable  accuracy.  The  fol¬ 
lowing  analyses  deal  with  approximate  methods. 

Coordinate  transformation 

A  method  suggested  by  Lin  (1971 )  uses  a  coordinate 
transformation  to  reduce  a  problem  with  a  variable  phase 
change  area,  such  as  a  cylinder,  to  one  with  a  constant 


Figure  3.  Accuracy  of  heat  balance  integral  solution,  phase  change  ST  =0.1,  a,  2  = 
k12  =  1.0. 
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phase  change  area-the  semi-infinite  solid.  Since  the 
relations  for  the  semi-infinite  solid  are  well  known  (see 
eq  1  and  2),  this  is  a  useful  procedure  but,  as  will  be 
seen,  is  limited  in  accuracy. 

The  following  transformation  will  reduce  the  cylind¬ 
rical  system  to  the  constant  area  case. 

y=r0\n(r/r0).  (34) 

The  phase  change  interface,  which  is  the  va'ue  of  y 
when  r  -  R,  is  related  by 

r)  =  r0  \n(R/r0).  (35) 

The  governing  equations  for  the  cylindrical  system 
then  transform  into  the  following  system  which  is 
valid  near  the  phase  change  interface  where  r  R. 

c£l  =  3  T(n,  r?)  3 T 

dy2  pta  3  y  3tj 

T(y,  0)  =  Tf 

T{v,v)  =  Tf 

r(o,7j)  =  rs. 

The  solutions  of  this  system  of  equations  are  universal 
functions  for  all  cross-sectional  areas.  However,  the 
solutions  are  only  valid  near  the  phase  change  interface. 
The  system  of  eq  36  need  not  be  solved  to  us*  the 
method.  The  phase  change  interface  rate  of  movement 
is  given  by 

dR  =A  \A^]2  9T(t 7,  v)  (37) 

dt  pi  [/((«)]  3 y 

For  the  constant  area  case  A  (r0)  =  A  (R)  and 

dJL=- 1  VkLiL  =g(n).  (38) 

dt  pi  3  y 

Thus  if  the  velocity  of  the  phase  change  interface  for 
the  constant  area  case  is  given  by 

d4  =  9W  (40) 

dt 

then,  the  phase  change  interface  velocity  for  the  cylin¬ 
drical  system  is 


(36) 

(36a) 

(36b) 

(36c) 


The  plane  or  Neumann  solution  is  given  by  eq  1  and 
2.  From  these  relations 


(42) 


Finally,  the  velocity  of  the  cylindrical  interface  is 


dR  -  2>2fti 
dt  R  \n(R/r0) 


(43) 


Equation  43  may  be  integrated  to  give 


2/J2  ln/J-02+l  =  *2-t. 

% 


(44) 


Equation  44  may  be  compared  to  the  quasi-steady  solu¬ 
tion,  eq  8,  for  the  case  of  no  superheat,  i.e.  ST  =  0, 

0  =  0.  In  the  limit  as  0-+O,  eq  2  is  -y2  =  (ST/2+5T) 
and  in  this  case  eq  44  is  identical  to  eq  8. 

While  eq  44  is  in  an  extremely  convenient  form  for 
cylindrical  systems,  its  accuracy  is  limited  to  certain 
ranges  of  ST,  0  and  t.  This  can  be  seen  by  comparing 
the  phase  change  interface  position  /3  for  a  cylindrical 
system,  calculated  with  eq  44,  and  the  values  computed 
numerically  by  Sparrow  et  al.  (1978),  given  in  Figures 
23  and  24. 

Equation  44  is  accurate  if  r/ST  <  1.0,  when  0  =  4. 
For  smaller  values  of  0  the  time  limit  when  eq  44  is 
accurate  will  increase.  Equation  44  is  so  simple  that 
it  may  be  of  value  for  quick,  more  or  less  crude,  esti¬ 
mates,  especially  when  0  0.0.  However  a  more  accu¬ 

rate  closed  form  relation  for  cylinders  with  superheat 
or  subcooling  will  be  evaluated  in  the  next  section. 


Effective  thermal  diffusivity 

Churchill  and  Gupta  (1977)  have  introduced  another 
method  which  allows  the  Neumann  solution  to  De  used 
for  more  complex  geometries.  The  procedure  involves 
replacing  the  nonlinear  phase  change  problem  with  its 
linear  analogue  which  does  not  include  phase  change. 
The  thermal  diffusivity  of  the  latter  problem  is  then 
replaced  by  an  effective  diffusivity  which  includes  the 
latent  heat.  Since  many  solutions  are  available  for 
non-phase  change  problems  the  method  has  potential 
for  application  to  numerous  freeze/thaw  problems. 

The  derivation  of  the  effective  diffusivity  to  use  is 
based  upon  the  fact  that  the  solution  to  the  zero  latent 
heat  analogue  of  the  Neumann  problem  (simply  transi¬ 
ent  conduction  in  the  semi-infinite  medium)  can  be 
forced  to  agree  with  the  Neumann  solution  if  an  effec¬ 
tive  diffusivity  replaces  the  actual  diffusivity.  If  the 
location  of  the  freezing  surface  is  important  then  one 
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change  position  vs  time  xe  =  0. 1.  F'Sure  7-  chan^  position  vs  time  xt  =  (X2. 


Figure  JO.  Phase  change  position  vs  time  xe  =  0.5.  Figure  1 1 .  Phase  change  position 


Freezt 


position  vs  time  xc  -  0.9,  freeze.  Figure  19.  Phase  change  position  vs  time  =  1.0  (pure  water),  thaw 


Figure  23.  Phase  change  position  vs  time,  Sj  =  0.10;  coordinate  transformation 
method. 


method. 

value  of  effective  diffusivity  is  used,  while  if  the  sur¬ 
face  heat  flux  density  is  desired  a  second  thermal  dif¬ 
fusivity  is  used. 

The  following  relations  can  be  used  for  the  heat  flux 
density  and  the  phase  change  location,  respectively. 

(45) 

(46) 
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where 


a,  = 


“thawed  for  a  thawing  problem 


“frozen  for  a  freezing  problem 


and  erf  is  the  error  function. 

The  method  consists  of  first  solving  the  linear  con¬ 
duction  problem  and  then  finding  the  relation  between 
the  space  coordinate  and  the  freezing  (or  simply  the 


32°F)  isotherm.  In  this  relation  one  replaces  the  ther¬ 
mal  diffusivity  by  eq  46  to  obtain  the  phase  change  in¬ 
terface  location  for  the  phase  change  problem.  It  is 
likely  that  numerical  evaluation  will  be  necessary  be¬ 
cause  of  the  complex  inverse  relation  between  the 
space  coordinate  and  the  freezing  isotherm  in  the  linear 
problem.  The  method  is  powerful  enough  to  handle 
complicated  problems  but  the  range  of  validity  is  not 
known.  The  method  cannot  take  into  account  the  dif¬ 
ference  in  thermal  properties  of  the  regions  for  which 
the  temperature  is  above  or  below  32°F  when  finding 
the  no-phase  change  solution. 

Churchill  and  Gupta  (1977)  applied  the  method  to 
cylinders  and  to  freezing  in  a  corner  with  good  results. 
They  used  the  exact  solution  for  the  cylindrical,  no¬ 
phase  change  problem  given  by  Carslaw  and  Jaeger 
(1959).  This  required  the  use  of  tabulated,  numerical 
values.  The  following  section  will  derive  a  convenient 
approximate  solution  to  the  problem  which  yields 
good  results. 

The  first  step  in  the  solution  of  the  phase  change 
problem  is  to  solve  for  the  temperature  of  the  pure 
conduction  problem  with  zero  latent  heat.  For  the 
case  of  a  cylinder  surrounded  by  an  infinite  medium, 
an  exact  solution  for  the  surface  step  change  situation 
is  given  by  Carslaw  and  Jaeger  (1959).  Flowever,  the 
solution  involves  an  infinite  series  of  Bessel  functions 
which  were  approximated  by  error  functions  for  small 
values  of  time.  The  final  results  are  in  graphical  form. 
An  approximate  solution  to  this  problem  can  be  found 
by  using  the  heat  balance  integral  methods  and  the 
equations  are  given  below,  referring  again  to  Figure  1. 
The  melting  case  will  be  examined  but  the  results 
apply  to  freezing  also. 


T 16,  t)  =  T0  (47a) 

T(r0,  t)  =  Ts  (47b) 

{S,  t)  =  0.  (47c) 

dr 


The  heat  balance  integral  of  eq  47  is  a  single  integration 
over  space  and  reduces  to 


Mrp,  t)  _  ^ 
dr  dt 


(48) 


0  = 


(49) 


where  v=  ( T~T0)/(TS-T0 ).  For  simplicity,  a  polyno¬ 
mial  in  r  is  assumed  for  the  temperature 

•'  =  (ff^)n:  «  (50) 

which  satisfies  eq  47a-47c  and  also  the  smoothness 
relations 

(51) 

dr"-1 

Boundary  layer  theory  has  shown  that  the  additional 
smoothness  relations  of  eq  51  may  improve  the  accu¬ 
racy  of  integral  methods  to  a  certain  extent.  The  ini¬ 
tial  condition,  v[r,  0)  =  0,  cannot  be  met  but  this  will 
not  seriously  affect  the  final  results. 

Equations  40  through  50  lead  to  a  differential  equa¬ 
tion  for  6  as  follows 

/2A2  =  "kill  (52) 

\/)  +  2  /  dr  ST 

where  A  =  (6-r0)/r0.  Equation  52  is  easily  integrated 
to  give 


—  4.  AJ+A2  =  2n(n  +  1)J_.  (53) 

3(n+2)  ST 


Volkov  and  Li-Orlov  (1 970)  noted  that  the  accuracy  of 
the  integral  method  could  be  improved  by  integrating 
eq  42  twice  over  the  space  coordinate.  El-Genk  and 
Cronenberg  (1979)  applied  this  idea  to  the  Neumann 
problem  with  apparent  success.  However,  for  the 
cylindrical  system  this  resulted  in  considerably  poorer 
results  than  the  integral  heat  balance  for  any  given  value 
of  n  >  2. 

Equations  50  and  53  complete  the  solution  of  the 
no-phase  change  problem.  These  equations  are  accept¬ 
able  when  compared  to  the  exact  results  of  Carslaw  and 
Jaeger  (1 959).  The  location  of  the  phase  change  inter¬ 
face  is  found  by  using  the  movement  of  the  isotherm 
with  the  fusion  value.  Thus  eq  50  gives 


v 


f 


(54) 


The  location  of  the  fusion  value  isotherm  is  then 

‘■■t&rj  (5si 

where  Aj  =  [rf-rj/r0.  Finally,  the  effective  thermal 
diffusivity  (eq  46)  will  replace  the  thermal  diffusivity 
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Figure  25.  Phase  change  vs  time,  accuracy  ST  =  0.10,  a12  -  k)2  =  I.Oof  ef¬ 
fective  diffusivity  method. 


in  eq  53.  The  solution  for  the  actual  phase  change  in 
terface  location  is  then 


_ 1 — 0*3  =  /3*2 -2/7(/?+1)r*  =  0 

3(n+2) 


where 


T* 


«eft_  act 

r2  O 


T 


(56) 


0+1 

The  accuracy  of  the  solution  increases  as  n  increases. 
Above  n  =  20  the  improvement  is  slight  and  thus  zt  =  20 
was  used  for  the  numerical  evaluations.  With  n  =  20  the 
solution  of  eq  56,  explicitly  for  0  as  a  function  of  time, 
is 


Equation  57,  combined  with  eq  46  for  the  effective 
thermal  diffusivity  and  eq  2  for  y,  is  a  simple  relation 
for  freezing  or  thawing  about  a  cylinder  which  is  accu¬ 
rate  enough  for  most  engineering  purposes.  A  compari¬ 
son  of  eq  57  with  previous  numerical  results  is  shown 
in  Figures  25  and  26.  These  figures  are  limited  to 
values  offc12  =(*i2  =  I.Obuteq  57  can  be  used  for 
any  ratios  of  the  frozen  and  thawed  properties.  It  is 
thus  considerably  more  flexible  than  the  numerical 
results.  For  most  cases  it  will  be  accurate  to  within 
10%  of  the  numerical  results.  This  is  felt  to  be  satis¬ 
factory  for  engineering  calculations. 


Heat  transfer 

The  instantaneous  heat  flow  from  the  cylinder  and 
the  total  heat  loss  or  gain  during  a  given  time  are  also 
quantities  of  interest.  The  surface  heat  transfer  rate  is 
given  by 


q  =  -61  A 


*nrvt) 

dr 


(58) 


P  =  [l  '(^7)°  05]f(<7+^1/3  +  (o-cf)1/3  -5.5]  +1 

(57) 


<7  =  6930-*?  -I- -166.375 
a  ST 

d=(a2- 27680.6406) 1/2 . 


Use  of  eq  50  and  57  leads  te  the  following  nondimen- 
sional  heat  flow 


2nk,S  e 


where  e  =  (o+d)1^3  +  {a-d)''^  -  5.5  as  before,  and  eq 
45  (aeq/cti )  =  [  (0+ 1 )  erfy) 2  is  used  with  a  and  d. 


) 
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Figure  26.  Phase  change  vs  time,  accuracy  ST  »  1.0,  <11,2  =  kl2  =  1-0  of  effective 
diffusivity  method. 


Figure  27.  Heat  flux  vs  time,  accuracy  ofa12  =  k]2  =  1-0,  Sy=0.1,  effective  diffusiv¬ 
ity  method. 


The  integrated  heat  transfer  at  the  cylinder  surface  Q ,  _  +A.  (61 ) 

is  ‘  2  nrip,t  21  '22  ) 


This  can  be  written  in  nondimensional  form  as 


In  eq  59  and  61,  the  value  n  =  20  has  again  been  used. 
Figures  27  and  28  show  q*  and  Q*  plotted  for  ST  = 

0.1  and  1 .0  with  <*12  =  *12  =  1  The  results  comPare 
fairly  well  with  published  numerical  values.  The  effec¬ 
tive  diffusivity  for  heat  transfer  does  not  seem  as  accu¬ 
rate  as  that  for  phase  change  location  but  it  still  gives 


16 


T 


T 


T 


T 


T 


ST  =  IO 

(•)  Tien  and  Churchill  C65) 
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.Eg  59  or  61 


Figure  28.  Heat  flux  vs  time,  accuracy  of  a12  =  k ■, 2  -  1.0,  ST  =  / .0;  effective 
diffusivity  method. 


reasonable  results.  Equations  59  and  61  can  be  applied 
to  any  values  of  a]2  and  &]2.  The  method  of  this  sec¬ 
tion  cannot  be  applied  when  the  superheat  parameter 
0  is  zero  but  for  this  case  there  are  sufficiently  accurate 
solutions  available,  as  given  earlier. 


The  problem  of  evaluating  the  effects  of  freeze/thaw 
about  a  cylinder  surrounded  by  a  homogeneous,  infinite 
medium  has  been  completed.  The  solutions,  while  not 
exact  (except  for  the  straight  numerical  results),  are 
sufficiently  accurate  for  engineering  design. 


CONCLUSIONS 
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have  also  been  derived,  given  by  eq  57  and  59.  These 

i  be  used  for  any  combination  of  ST,  0,  <*]2,  612 
and  will  be  accurate  enough  for  most  permafrost  cal¬ 
culations.  These  solutions  are  limited  to  cases  for 
which  0  >  0.  If  the  superheat  is  zero,  then  eq  1 4  or 
eq  44  can  be  used  with  good  results.  The  effective 
thermal  diffusivity  method  used  here  should  be  adapt¬ 
able  to  many  problems  of  freezing  and  thawing. 
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